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Abstract. We report results from numerical simulations of star formation in the early universe 
that focus on the role of subsonic turbulence, and investigate whether it can induce fragmentation 
of the gas. We find that dense primordial gas is highly susceptible to fragmentation, even for rms 
turbulent velocity dispersions as low as 20% of the initial sound speed. The resulting fragments 
cover over two orders of magnitude in mass, ranging from ~ 0.1 M© to ~ 40 M©. However, 
our results suggest that the details of the fragmentation depend on the local properties of the 
turbulent velocity field and hence we expect considerable variations in the resulting stellar mass 
spectrum in different halos. 
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1. Introduction 

Over the course of the last decade, work by a number of groups has lead to the 
development of a widely accepted picture for the formation of the first stars, the so- 
called Population III (or Pop. Ill) stars. In this picture, the very first stars (often termed 
Population III.l) form within small dark matter halos that have total masses M ~ 
10^ M0, virial temp eratures of around a thousand K, and that are assem bled at redshifts 
z ~ 25-30 or above (JBromm fc Larsonl20o4lGlovei<2o6HBromm et al.|[2009il . Gas falhng 
into one of these small dark matter halos is shock-heated to a temperature close to the 
virial temperature of the halo, and thereafter cools via II2 rotational and vibrational 
line emission. A cold, dense core forms at the center of the dark matter halo, with a 



characteristic density ?ich ^ 10 cm^'' and temperature Teh ^ 200 K; these values are 
determined by the microphysics of the H2 molecule, and represent the point at which 
II2 cooling becomes inefficient, and the effective equation of state of the gas becomes 
stiffer than isothermal. Once the mass of this cold dense clump exceeds the critical 
Bonnor-Ebert mass for this temperature and density, which is approximately 1000 Mq, 
it decouples from the larger scale flow and undergoes runaway gravitational collapse. Until 
rece ntly, numerical simulations that have followed this coll apse up to very high densities 
(e.g. lAbel. Bryan, fc NormanI l2002t iBromm fc Loebl 12004 1 have found no indication of 
further fragmentation, and so it has typically been assumed that this dense core forms 
only a single star. In this picture, the final mass of this star will be set only by the 
efficiency with which it can accrete gas from its natal core, and there are good reasons 
to believe that for protostellar masses less than ^ 20-30 M0, this efficiency will be close 
to 100% ( McKee fc TanI 120081 ). It is therefore difficult to avoid the conclusion that the 
final mass of the Pop. Ill star will be very large. 

Over the past few years, however, a new generation of numerical simulations have be- 
gun to suggest that the true picture of Pop. Ill star formation may be rather more 
com plicated than the simp le scenario outlined abo ve. Work by a number of groups 
Clark. Glover fc Klessenl2008.: Machida 2008: Machida et al.l2Q09HTurk. Abel fc O'Shea 



2009tlStacv. Greif fc Brommll2010 ) has found evidence for fragmentation in high density 
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primordial ga s, with predictions of the like lihood that fragment ation occurs ranging from 
roughly 20% (iTurk. Abel fc O'Sheal l20"Q9h to essentially 100% (IClark. Glover fc Klessen 
2008tlStacv. Greif fc Bromnj|2010l) . If fragmentation, and the formation of binary or mul- 
tiple systems, is indeed a common outcome of Population III star formation, then this 
has profound implications for the final masses of the Pop. Ill stars, their production of 
ionizing photons and metals, the rate of high redshift gamma ray bursts, and many other 
issues. 

It is therefore important to better understand the physical basis of fragmentation in 
these systems. Although numerical simulations starting from cosmological initial condi- 
tions and following collapse all the way to protostellar densities will give us the most 
accurate picture of how the first stars formed, these simulations are computationally 
costly (making it difficult to explore a wide parameter space) , and are also difficult to in- 
terpret, since many physical effects combine to yield the observed behaviour. Simulations 
that start from simpler initial conditions and that focus on exploring the importance of 
a single free parameter therefore play an important role, which is complementary to that 
of fully cosmo logical simulat i ons. A good exampl e is the recent work by Machida and 
collaborators ( Machidall2008t iMachida et al.ll2009l ) which examined the influence of the 
initial rotational energy of the gas, and found that zero metallicity clouds with sufficient 
initial rotational energy could fragment into tight binaries. In this contribution, we briefly 
describe the results of a similar study that we hav e recently perform ed that looked at 
the effects of the initial turbulent energy of the gas ( Clark et al.ll2010l) . 



2. Simulations 



2.1. Numerical method 



We modeled the evolution of the primordial gas using a modified version of the Gadget 
2 smoothed particle hydrodynamics code (Springcl 2005). Our modifications in clude the 



addition of a sink particle implementation ([Bate. Bonnell fc Pried Il995t Ijappsen et al. 



20051) ■ allowing us to follow the evolution of the gas beyond the point at which the 
first protostar forms , and a detailed treatment of primordial coo hng and ch e mistr y 
(JGlover fc AbellbOOsHcTark et al.ll2010l ). FuU details can be found in lClark et all (|2010l ). 



2.2. Initial conditions 

We took initial configurations that were unstable Bonnor-Ebert spheres, into which we 
injected a subsonic turbulent velocity field. In our study, we considered two different sets 
of initial conditions: one corresponding to the Pop. III.l scenario described above, and 
another corresponding to the Pop. III. 2 scenario, in which efficient HD cooling leads to a 
lower Teh and higher rich for the primordial prestellar core. However, in this contribution 
we will discuss only the results of the Pop. III.l simulations. In these, we set the initial 
gas temperature to be 300 K, and took an initial central number density of 10^ cm~^. 
The mass of the Bonnor-Ebert sphere was taken to be IOOOM0, corresponding to roughly 
three Jeans masses. This mass was represented by 2 million SPH particles, giving us a 
particle mass of 5 x lO^'' Mq, and a mass resolution of 5 x 10^^ M©. Within the BE 
spheres, we imposed a turbulent velocity field that had a power spectrum P{k) oc k~ . 
The three-dimensional rms velocity in the turbulent field, Auturb, was scaled to some 
fraction of the initial sound speed Cg . For the simulations presented here, we used four 
different rms velocities: 0.1, 0.2, 0.4 and 0.8 Cg. Sink particles were created once the 



number density of the gas reached 10 cm ^, using the standard iBate. Bonnell fc Price 



( 19951) algorithm, with a sink accretion radius of 20 AU 
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Figure 1. The top panel shows the mass functions from the simulations in which fragmentation 
occurs. In all cases the mass function is plotted at the point at which the total mass of gas 
converted to sink particles is lOOM©. Note that as accretion is ongoing, and the system is still 
young (t ~ 1000 yr), these will often not be the final masses of the sinks. The mass functions in 
the individual simulations differ substantially, although the combined mass function, shown in 
the bottom panel, exhibits a broad and fiat distribution between masses of 0.4 and 20 M©. 

3. Key results 

We find that the clouds with A?;turb ^ 0.2cs all fragment into small clusters of sink 
particles: the runs with Awturb = 0.2cs, 0.4cs and 0.8cs have formed 5, 31, and 15 sink 
particles, respectively, at the point at which 10% of the cloud mass has been accreted. This 
clearly demonstrates that turbulent subsonic motions are able to promote fragmentation 
in primordial gas clouds. Only in the case with Auturb = O.lcg does the cloud form 
just a single sink particle. In this case, an extended disk builds up around the star, but 
otherwise the gas contains little structure. 

Interestingly, there is no clear trend linking the number of fragments that form and 
the initial turbulent energy: the 0.4 Cg run fragments more than the 0.8 Cg run, despite 
containing only one quarter of the initial turbulent energy. The effects at play here are 
somewhat complex. First, the turbulent velocity field contained in the collapsing region 
will differ with each value of Awturb/cs, since the different strengths of flow will push the 
gas around to different degrees. In addition, the nature of the turbulence that survives 
in the collapsing core will also affect the fragmentation. The cloud with Auturb = 0.4cs - 
the most successful in terms of fragmentation - forms a large disk-like structure, which 
appears to aid the fragmentation of the infalling envelope. Much of the cloud's ability to 
fragment therefore seems to depend on the amount of angular momentum that becomes 
locked up in the collapsing region, and hence on the detailed properties of the turbulent 
velocity field. Further, the extra delay in the collapse caused by the increased turbulent 
support also gives the cloud more time to wash out anisotropics. The ability of the 
cloud to fragment is a competition between these conflicting processes. For the randomly 
generated velocity field that we used in this study, the ability to fragment is better for 
Awturb/cs — 0.4, than Auturb/cs — 0.8, but we note that this may not always be the case. 

The mass functions of the sink particles from the simulations that undergo fragmenta- 
tion are shown in Fig.[TJ For clarity, we have omitted the single 100 M© sink particle that 
forms in the 0.1 Cg cloud. For the 0.2 Cg and 0.8 Cg clouds, we see that the sink masses 
cluster around some central value - roughly 12 M© and 4 Mq respectively - while the 0.4 
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Cs cloud has a mass function that is skewed to lower masses, with a peak at around 1 Mq, 
a sharp fall-off below this, and a broad distribution towards higher masses, extending 
up to around 13 Mq. As all of the sinks form with essentially the same initial mass, 
the spread of sink masses at the end of the simulation reflects the fact that some sinks 
are more effective than others at accreting gas from the available reservoir within the 
cloud. This behaviour is very familiar from studies of present-day star formation (see e.g. 
Bonnell. Vine fc Batell2004tlSchmeia fc Klessenll2004 ). and is a result of the the high sen- 



sitivity of the mass accretion rate to the sink mass and the relative velocity between sink 
and gas (JBonnell et al.ll2001al ): r/i* oc pm'l/v'^^Y, where to* is the mass of the sink particle. 



p is the gas density, and Urei is the velocity of the sink particle relative to the gas. As the 
number of sinks increases, more and more get 'kicked' by close dynamical interactions, 
forcing them onto more distant orbits, and increasing Wrci. These sinks therefore accrete 
very little gas following the dynamical kick. The few sinks that are not strongly kicked 
therefore accrete the lion's share of the available gas, and form the high- mass tail of the 



resulti ng mass function. This process is termed 'competitive accretion', (jBonnell et al. 



2001al) and normally leads to the type of distribution of masses seen in our 0.4cs run. 

Our simulations therefore demonstrate that if significant levels of subsonic turbulence 
are present in the primordial gas found within the first dark matter halos, then Pop. Ill 
stars forming in these halos are more likely to be born in small stellar groups than as 
single stars. Moreover, if Pop. Ill stars do indeed form in this way, then we would expect 
the resulting initial mass function to be broad, thanks to the influence of competitive 
accretion within the star cluster. Tests of these expectations will require large, high- 
resolution studies of the formation of the first dark matter minihalos that are capable 
of resolving the turbulent flows in the gas at the moment when the baryons become 
self-gravitating, and that use sink particles (or some comparable technique) to allow the 
evolution of the self-gravitating system to be followed over multiple dynamical times. 
Such simulations would help determine which, if any, of the initial conditions presented 
in our study are actually realized in nature. 
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